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Abstract 

We have performed self-consistent 2.5-dimensional nonsteady MHD numerical sim- 
ulations of jet formation as long as possible, including the dynamics of accretion 
disks. Although the previous nonsteady MHD simulations for astrophysical jets re- 
vealed that the characteristics of nonsteady jets are similar to those of steady jets, 
the calculation time of these simulations is very short compared with the time scale of 
observed jets. Thus we have investigated long term evolutions of mass accretion rate, 
mass outflow rate, jet velocity, and various energy flux. We found that the ejection 
of jet is quasi-periodic. The period of the ejection is related to the time needed for 
the initial magnetic filed to be twisted to generate toroidal filed 

Tejection j— OC — OC Emg ■ 
V A D 

We compare our results with both the steady state theory and previous 2.5- 
dimensional nonsteady MHD simulations. Then it is found that time averaged velocity 
of jets (Vj e t ave) are ~ O.IVk and ~ 0.1Vj et max , where Vk is the Keplerian velocity at 
(r,z) = (1,0) and Vj etiaax is the maximum velocity of jet. Nevertheless, the charac- 
teristics of our simulations are consistent with those of steady solution and previous 
short term simulations in that the dependences of the time averaged velocity V z ave 
and mass outflow rate M w ave on the initial magnetic field strength are approximately 

M w ave ex B - 32 and V jet ave ex .Bo ' 3 . 
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1. Introduction 

Astrophysical jets have been observed in young stellar objects (YSOs) (e.g., Fukui et 
al 1993; Ogura 1995; Burrows 1996; Bachiller 1996; Reipurth et al. 2002; Curiel et al. 
2006), active galactic nuclei (AGNs) (e.g., Bridle & Perley 1984; Biretta et al. 1995; Junor 
et al. 1999; Jiang & Hong 2003; Matthew 2005), and some X-ray binaries (XRBs) (e.g., 
Margon 1984; Mirabel et al. 1994; Tingay et al. 1995; Kotani et al. 1996; Migliari et al. 
2006). Although the acceleration and collimation mechanisms of these jets are still not well 
understood, these objects are believed to have accretion disks in their central regions. 

The standard model of a jet-disk system has been an original idea of Blandford & Payne 
(1982). Energy and angular momentum are removed magnetically from the accretion disk by 
field lines anchored to the disk surface and extending to large distances. Blandford & Payne 
(1982) showed that a centrifugally driven outflow of matter from the disk is possible, if the angle 
between field line and disk is less than 60°. They discussed self-similar solutions of the steady 
and axisymmetric MHD equations and the possibility of such acceleration and collimation of 
the flow from a cold Keplerian disk. After their work, many authors studied MHD models of 
jet formation from accretion disks (e.g., Pudritz & Norman 1986; Sakurai 1987; Lovelace et al. 
1991; Contopoulos & Lovelace 1994; Najita & Shu 1994; Fendt & Camenzind 1996) based 
on the theory of steady and axisymmetric MHD winds, which was first developed by Weber 
& Davis (1967) for the solar wind. Cao & Spruit (1994) examined the mass outflow rate of 
magnetically driven jets by studying the solution that passes through the slow magnetosonic 
point (see also Li 1995). They confirmed that the inclination angle of the field line is very 
important for achieving a high mass outflow rate. 

Kudoh & Shibata (1995, 1997a) studied one-dimensional steady magnetically driven jets 
along a fixed poloidal field line for a wide range of parameters, assuming the shape of that field 
line. They found that the steady solutions generally can be classified into two branches: (1) 
magneto-centrifugally driven jets, when the magnetic field is strong, in that case they found 
that the jet velocity (Vj et ) and mass outflow rate (M w = ^jf 1 ) depend on the magnetic energy 
(E mg ) as Vj e t oci?o 1//3 , M w ~ constant and (2) magnetic pressure-driven, when the magnetic field 
is weak, in that case Vj e t oc So 176 , M w oc .Bo ' 5 . In the centrifugal case, the main acceleration 
occurs before the flow speed exceeds the Alfven speed, while in the latter case most of the 
acceleration occurs after the flow speed becomes greater than the Alfven speed. This mass flux 
scaling law was confirmed by Ustyugova et al. (1999). 

In most theoretical models of jets from accretion disks, however, accretion disks are 
treated as boundary conditions. Accretion disks only play a role of supplying energy and 
mass to the jets, and neither accretion flow nor internal structure of disks are considered in 
these models. Since the disk itself is not treated, such simulations may last over hundreds of 
Keplerian periods. This idea was first applied by Ustyugova et al. (1995). Extending this work, 
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Romanova et al. (1997) found a stationary final state of a slowly collimating disk wind in the 
case of a spilt-monopole initial field structure after 100 Keplerian periods. Ouyed & Pudritz 
(1997a, 1997b) presented time-dependent simulations of the jet formation from a Keplerian 
disk. For a certain (already collimating) initial magnetic field distribution, a stationary state of 
the jet flow was obtained after about 400 Keplerian periods of the inner disk with an increased 
degree of collimation. Ouyed et al. (2003) investigate the problem of jet stability and magnetic 
collimation extending the axisymmetric simulations to 3D. 

On the other hand, Uchida & Shibata (1985) and Shibata & Uchida (1987, 1989, 1990) 
performed time-dependent, two-dimensional axisymmetric (2.5-dimensional) MHD numerical 
simulations of magnetically driven jets from accretion disks. They solved the interaction be- 
tween a geometrically thin rotating disk, including the dynamics of the disk, and a large-scale 
magnetic field that was initially uniform and vertical . Shibata & Uchida (1986) investigated 
the detailed properties of these jets. They found that (1) the velocity of the jet was typically 
of order of the disk's Keplerian velocity and (2) it increased with increasing magnetic field 
strength in a manner similar to the scaling law of Michel's (1969) solution. Matsumoto et al. 
(1996) carried out 2D MHD simulations of a torus threaded by poloidal magnetic fields and 
found that (1) the jet velocities were again comparable to the Keplerian velocity and (2) the 
mass outflow rate of the jet increased with the strength of the initial magnetic field. Kudoh 
et al (1998) studied the formation mechanism of jets from geometrically thick disks and the 
dependence of the initial magnetic field strength (B ) in detail by performing self-consistent 
2.5-dimensional , nonsteady, ideal MHD simulations including the dynamics of the disk. They 
found that the ejection point is determined by the effective potential resulting from the gravi- 
tational and centrifugal forces along the field lines (Blandford & Payne 1982) and also that the 
velocity and mass outflow rate are consistent with those predicted by the steady theory (Kudoh 
& Shibata 1995 and Kudoh & Shibata 1997a). Kato et al (2002) performed 2.5-dimensional, 
axisymmetric, ideal MHD simulations of jets from geometrically thin disks for Keplerian and 
sub-Keplerian cases over a wide range of initial magnetic field strengths. Kigure & Shibata 
(2005) investigate the problem of jet formation and stability by using 3-dimensional MHD sim- 
ulations. To investigate the stability of the MHD jet, they introduce a perturbation to the 
accretion disk with a nonaxisymmetric sinusoidal or random fluctuation of rotational velocity. 
In both perturbation nonaxisymmetric structure with m = 2 appears in the jet, where 

m is the azimuthal wavenumber. They conclude that this structure seems to originate in the 
accretion disk. 

On the other hand, the acceleration and collimation of the jet have been studied in 
the steady (e.g. Sauty et al. 2004; Bogovalov & Tsinganos 2005). In recent years, many 
simulations taking other physical processes into consideration, e.g. the magnetic diffusion ( 
Kuwabara et al. 2000, Kuwabara et al. 2005; Fendt & Cemeljic 2002; Casse & Keppens 
2002, Casse & Keppens 2004), the dynamo process in the accretion disk (von Rekowski et al. 
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2003), and the radiation force ( Proga 2003), have been performed. Koide et al. (1999) showed 
that a magnetically driven jet in the general relativistic MHD simulation has characteristics 
similar to those of the nonrelativistic MHD jet (Shibata & Uchida 1986). In addition to these 
studies which consider the initial magnetic field as large scale, several papers were considering 
the evolution of a stellar magnetic dipole in interaction with a diffusive accretion disk. Hayashi 
et al. (1996) observed magnetic reconnection and evolution of X-ray flares during the first 
rotational periods in their numerical simulations. 

We note that all calculations including the treatment of the disk structure and the ideal 
MHD, (e.g. Shibata & Uchida 1986, Shibata & Uchida 1987, Shibata et al. 1989, Shibata 
& Uchida 1990; Matsumoto et al. 1996; Kudoh et al. 1998; Kato et al. 2002; Kudoh et 
al. 2002 and Kigure & Shibata 2005), could be performed only for 1 — 2 Keplerian periods of 
inner disk. At this point, we emphasize that the observed kinematic time scale of protostellar 
jets can be as large as 10 3 — 10 4 yrs, corresponding to 5 x 10 4 — 5 x 10 5 stellar rotational periods 
(and inner disk rotations)! For example, proper motion measurements for HH30 jet (Burrows 
1996) give a knot velocity of about 100 — 300/cms -1 and a knot production rate of about 0.4 
knot per year. Assuming a similar jet velocity a long the whole jet extending a long 0.25 pc 
(Lopez et al. 1995), the kinematic age is about 1000 yrs 

In order to have access to the observed time scale of jets and to know whether the jet for- 
mation becomes quasi-steady state, we performed long term 2.5-dimensional MHD simulations 
of jets. 

We solve the dynamics of the disk as in Shibata & Uchida (1986), Kudoh et al. (1998) 
and their following studies. We also want to know whether the time averaged physical quantities 
have the same characteristics as those in the steady model and previous simulations (Kudoh et 
al. 1998) and (Kato et al. 2002). Our calculation time is about 20 times longer than those 
of previous simulations. Also, the simulation box is large enough to minimize the effect of top 
and side boundary conditions 

2. Numerical Method 

2.1. Assumptions 

Our simulations make the following assumptions: (1) axial symmetry around the rota- 
tional axis (d/dip = 0), including the azimuthal components of a velocity (v v ) and a magnetic 
field {B^) (i.e., 2.5-dimensional approximation); (2) ideal MHD (i.e., Ohmic diffusivity and/or 
ambipolar diffusion are assumed to be negligible); (3) an inviscid perfect gas with a specific heat 
ratio of 7 = 5/3; (4) a point-mass gravitational potential only, with disk self-gravity neglected. 

2.2. Basic Equations 

The basic equations we use are the 2.5-dimensional ideal MHD equations in cylindrical 
coordinates (r,ip,z) : 
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where \& = —GM / (r 2 + z 2 ) 1 / 2 is the gravitational potential, G is the gravitational constant, and 
M is the mass of a central object. Other variables are summarized in Table 1. 

2.3. Initial Conditions 

As an initial condition, we assume an equilibrium disk rotating in a central point-mass 
gravitational potential (Matsumoto et al. 1996). Exact solutions for these conditions can 
be obtained under the simplifying assumptions for the distribution of angular momentum and 
pressure (e.g., Abramowicz et al. 1978): 

L = L r a , (10) 

p = Kp l+1 ' n . (11) 

Then the distribution of material in the disk is given by 
GM 1 2 2 „ 2 p 

P 



constant. 



(r 2 + z 2 )V2 1 2(1 - a) ~ u ' v " ' """" ^ 

We use a = 0.45 and n = 3 throughout this paper. These are the different parameters from 
the previous similar simulations (e.g.,Kudoh et al. 1998, Kato et al. 2002). Our parameters 
make the disk thicker and wider than that of Kato et al. (2002), mainly for both numerical 
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p = p c exp< a 



(13) 



convenience and long term simulation. The mass distribution outside the disk is assumed to be 
that of a uniformly high temperature corona in hydrostatic equilibrium without rotation. The 
density distribution in hydrostatic equilibrium is 

r _ 

where r is a radius defined later (see Table 1), a = (jVko 2 /V SC 2 ),V SC is the sound velocity in the 
corona, Vko = (GM/ro) 1 / 2 is the Keplerian velocity at radius r . We use a = 1 and p c / Po = 10 -3 
throughout this paper, where p is a density defined later (see Tablel). 

For simplicity, the initial magnetic field is assumed to be uniform and parallel to the 
axis of rotation 

B Z = B = constant, B r = B^ = 0. (14) 

Indeed, the presence of a large-scale magnetic field in accretion- disk-jet systems is ob- 
served in AGN jets and protosteller jets (Ray et al. 1997; Pushkarev et al. 2005; Gabuzda 
et al. 2004; Wouter et al. 2006) and is assumed theoretically in jet launching models (e.g. 
Casse & Keppens 2004, and references therein) 

2.4- Boundary Conditions 

On the axis (r = 0) we assume a boundary condition that is symmetric for p,p,v z , and 
B z while v r ,v v ,B r , and B^ are antisymmetric. 

The side, top and bottom surfaces are free boundaries. In order to avoid a singularity 
at the origin, the region around r = z = is treated by softening the gravitational potential as 

-GM/{r 2 + z 2 f/ 2 for e < (r 2 + z 2 ) 1 / 2 , 

^ = | -GMl/e-[{r 2 + z 2 y/ 2 -e}/e 2 for 0.5e < (r 2 + z 2 ) 1 ' 2 < e, (15) 
^ -1.5GM/e for (r 2 + z 2 ) 2 < 0.5e. 

We use e = 0.2ro throughout this paper. 

These boundary conditions are different from those MHD simulations of jets from ac- 
cretion disks performed by Ouyed et al (1997a), Romanove et al. (1997), Meier et al. (1997), 
Ustyugova et al. (1999), and Pudritz et al.(2006). The main difference is the condition on the 
equatorial plane. They assumed an inflowing fixed boundary condition (i.e., the angular mo- 
mentum is continually injected from the boundary). Therefore, their numerical simulations did 
not allow disk accretion due to magnetic braking nor growth of magneto-rotational instability. 

2.5. Numerical Schemes 

The numerical schemes we use are the cubic-interpolated pseudo-particle (CIP) method 
(Yabe & Aoki 1991; Yabe et al. 1991) and the method of characteristics-constrained transport 
(MOCCT) (Evans & Hawley 1988; Stone & Norman 1992). The magnetic induction equation 
is solved using MOCCT and the others using CIP. This is the two-dimensional version of the 
scheme used by Kudoh & Shibata (1997b), and Koide et al. (1999) 
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2. 6. Parameter 

We normalize the physical quantities with their initial value at (r, z) = (tq, 0) taking 
r = (L 2 /GM) 1 " 1 " 2 "' so that the initial density of the disk has maximum value of po a t 
(r, z) = (ro,0). The normalized unit for each variable is summarized in Table 1. There are two 
nondimensional parameters: 

£th = ^4„ (16) 

^mg = — 2" (^) 
VKO 

where V s0 = (7Po/Po) 1//2 ; Kio — -Bo/(47rp ) 1//2 , and p is the initial pressure at (r, z) = (r ,0). 
The initial parameters we use in this paper are summarized in Table 2. We use E t h = 0.018 
throughout this paper. When E t ^ = 0.018 and a = 1, the temperature of corona is about 10 2 
times greater than that of the disk at (r, z) = (tq,0). In this paper, we adopt L = Lor OA5 
(Z/ = 1.00) that makes the disk thicker and wider than that of Kato et al. (2002), mainly for 
both numerical convenience and long term simulation: the internal structure of a thick disk 
is better resolved than that of a thin disk, and much matter in the disk can make the time 
of simulations longer than that of Kudoh et al. (1998). Of course, many accretion disks in 
YSOs, AGNs, and XRBs are expected to be geometrically thin and nearly Keplerian. Kato et 
al. (2002) also performed simulations for geometrically thin disk case. However they calculated 
only one or two orbits just the same as Kudoh et al. (1998). The minimum grid size is 0.05r in 
r-direction and O.Olro in z-direction. The maximum grid size in r-direction is O.lro and 0.4ro in 
z-direction, and the number of grid points used in this paper is 269 x 639. The grid spacing is 
uniform for r /ro < 1 and z/ro < 1 and stretched in r and z for r /r^ > 1 or z/zq > 1 . The size of the 
computational domain is (-R m i n /r ~-Rmax/ r 0; Z min /r ~ Z max /Vo) = (0.0 ~ 24.8, — 65.75~65.75). 

3. Numerical Results 

3.1. Typical Case 

Figures 1, 2 and 3 shows the time evolution of the density distribution (color scales, top 
panels), temperature distribution (color scales, bottom panels) (T = jp/p), poloidal magnetic 
field lines, and poloidal velocity (hereafter, the variables are expressed in nondimensional form) 
for model 1,4 and 6 respectively. Time t = 2tt ~ 6.28 corresponds to one Keplerian orbit at 
(r, z) = (1,0). Evolutions for all models look like similar each other. 

Now we show the case of model 4 (see Figure 2). In the early stage of evolution, a 
torsional Alfven wave is generated at the disk surface and propagates up into the corona (t = 13). 
Since this wave extracts angular momentum, the rotating disk begins to fall into the central 
region. The surface layer of the disk falls faster than the equatorial part. This is the avalanche- 
like accretion that was studied by Matsumoto et al. (1996). Because only a small fraction of the 
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accreted matter is ejected into the bipolar direction due to Lorentz force in the relaxing twist, 
both density and pressure of the inner region increase. Within this region, angular momentum 
is transferred from the high to the low-density parts. The cold material on the disk surface is 
ejected as a jet (t = 35). At t = 82 due to the magnetorotational instability the channel flow 
becomes clear. Initially both gas pressure and magneto-centrifugal force drives and accelerates 
the outflow (below the Alfven surface). After that, when the toroidal field generated as a result 
of the differential rotation of the accretion disk is accumulated, the acceleration is due to the 
magnetic pressure gradient, d(B^/87r)/dz. To illustrate the twist level of the magnetic field 
lines, we show in figure 4 the time variation of the ratio of the toroidal magnetic field, B v , to 
the poloidal magnetic field , B p , along r = 0.225. Figure 4 illustrates that in case of initial 
weak magnetic field the twisting field B^ become more dominant and both of the Alfven point 
and the slow point become near to the accretion disk (Pelletier & Pudritz 1992 and Kudoh & 
Shibata 1997a). In both the weak and strong initial magnetic field a strong twist appears at 
t=15 and propagates outward when the jet is ejected. 

The outflow consists of both the material that is initially in the corona and the material 
from the disk. Channel flow continues to grow in the disk and jets are ejected continuously and 
intermittently (t = 82). Jet ejection and accretion still continue at the last stage of evolution 
in our simulations (t = 115). We can also see that magnetic field lines entwine each other and 
that the magnetic turbulent flow develops. We can see that the magnetic islands are created by 
magnetic reconnections in the disk . Since we assumed ideal MHD, the magnetic reconnection 
is caused by numerical diffusion. Numerical diffusion is inevitable for finite difference numerical 
schemes, even though we do not include the magnetic diffusivity explicitly. The time evolution 
is similar to that in Matsumoto et al. (1996) and Kudoh et al. (1998). 

The avalanche-like flow is caused by a physical mechanism similar to that of the magne- 
torotational instability (Balbus & Hawley 1991; Hawley & Balbus 1991). Note that although 
the disk is initially rotating in our model, the corona is not. The discontinuity of azimuthal 
velocity at the interface between them generates torsional Alfven waves propagating into both 
the corona and the disk. Since the Alfven wave behaves like a large amplitude perturbation, it 
triggers the magnetorotational instability in the disk. The wavelength of this instability is deter- 
mined by the most unstable wavelength (A) for the magnetorotational instability, A ~ 2ttVa/Q, 
so the base height of the avalanche depends on the initial magnetic field strength. 

4. Time Evolutions of Some Physical Quantities of Jets 

One of the most important purpose of this study is to clarify whether the jet ejection 
becomes steady state. In figures 5,6 and 9 we show the time evolutions of the mass outflow rate 
(M w ), the mass accretion rate (M a ), and the toroidal magnetic field energy {E mgt ) for model 1 
{E mg = 2 x 1(T 4 ), model 4 (E mg = 2 x 10 -8 ) and model 6 [E mg = 2 x 1(T 6 ). M w is defined by 
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where z v is z-coordinates for jets to pass through. The mass accretion rate for various models. 



are very nonsteady, and sometimes the ejected matter falls back onto the disk, (2) when initial 
magnetic field is weak, the jet ejection point goes away from the equatorial plane after long 
time simulation. Kudoh et al. (2002) and Kuwabara et al. (2000) found that mass accretion 
and mass ejection take place intermittently in the case without perturbation in the disk. The 
simulation time of both of them ranged from one orbit in case of Kuwabara et al. (2000) and 
three orbits in case of Kudoh et al. (2002). Our simulations in this paper more than 20 orbits 
for different models. From Figure 5, it is clear that the mass ejection flux is still intermittent 
until the last stage of evolution. The intermittent ejection is clear in all models in the figure 
5 but there are some differences between them. The first one, in case of strong magnetic field 
the intermittency is almost similar at the beginning of the simulation to that at the end of 
simulation. But in case of weak field, as the simulation goes on, the ejection mass flux becomes 
very nonsteady. Not only the intermittency increases in case of weak magnetic field but also 
the absolute amount of mass flux increases. But in case of strong magnetic field we notice the 
opposite. Figure 6 shows the mass accretion rate for models 1, 4 and 6. The mass accretion, like 
mass ejection, is intermittent until the last stage of our simulations. In both weak and strong 
initial magnetic field the mass accretion is intermittent. The general trend of intermittency 
of the mass accretion is similar to the intermittency of mass ejection. That similarity is more 
clearer for example in model 4 than other models. In model 4 the absolute value of mass 
accretion at the beginning of simulation is smaller than that at the end. We notice the trend 
in the mass ejection, also the same trend happen in the other models. 

Because of the mass accretion rate is highly variable, plotting it at one radius in time 
shows very little of the overall characters of the accretion. We therefore show values that are 
averaged over space or time, in order to obtain a better understanding of the overall accretion 
within the disk. 

Figure 7 shows the accretion rate against time, averaged between r = 0.52 and r = 10. 
This region is chosen since the accretion has little value at larger radii. This shows that even 
though the accretion rate is highly intermittent, it is predominantly positive. The behavior of 
the accretion is more complex in case of weak initial magnetic field after 10 orbits, because the 
evolution of the disk take a long time until it becomes more turbulent. The accretion, in case 
of strong initial magnetic field, is positive during the first 5 orbits. During that time, part of 
the subtracted angular momentum is ejected as a jet. The other part is transported at large 



M a is defined by 




(19) 



It is not easy to define (M w ) and (M a ), because of following reasons: (1) ejection and accretion 
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distances by the magnetic stresses. This is exactly what we see in Figure 8 , which shows the 
accretion rate as a relation of radius of the disk, averaged over the first 10 orbits. Accordingly, 
we see that the accretion is positive in the inner part of the disk and the accretion is negative 
beyond the radius 2-4. The negative accretion continue about 5 orbits. The accretion becomes 
positive again because the continuous ejection of the jet which means continuous subtraction 
of the angular momentum. 

Figure 9 shows the toroidal magnetic energy for model 1, 4 and 6. The toroidal magnetic 
field energy E mgt for various models is defined by 

^=-(r^^ + jT%^)A/. ,ar *- « 2o) 

The intermittency is clear and the trend is also similar to mass ejection. Hence we try to find 
the relation between mass accretion, ejection and magnetic field. The ejections are intermittent 
and seem to have periods for all models, and the periodicity seems different in the models with 
different initial magnetic field. Next we check the times of ejection peaks for different models 
and study the relation between that peak times with the time needed for the toroidal field to be 
accumulated after that untwisted in vertical direction carrying the mass flux. From the initial 
conditions equation 14, B z = B = constant, B r = B v = 0. With the rotation of the disk by 
angular velocity Q, the toroidal field B^ is generated from B z . 

VLB Z . (21) 



dT 

By integration w.r.t. time T then 

B v « QB Z T. (22) 
The ejection occurs when the magnetic energy equals the rotational energy, then 

a ~ i pV 2 ~ i pV 2 (23) 

Combining Eqs. 22 and 23, yields 

B^T^ 1 r, 1 o , 

-^-«2^ W 2^« (24) 

11-1 

Tejection OC — OC E m g , (26) 

V A D 

where Va is the Alfven speed, B is the initial magnetic field strength and T e j ect i on is the time 
corresponding to the peak mass ejection rate as shown in Figure 10. Figure 11 shows the 
relation between the initial magnetic field strength and the average time interval between the 
peaks of the mass ejection rate of models 1, 3, 4 and 6. Here the squares correspond to our 
numerical values and the solid line is the best fit to them with the following dependence 



average,eje ^ ^mg ' 
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The dotted line corresponds to the analytical relation equation 26. We notice that the de- 
pendence of the ejection time on the initial magnetic field (E mg ) is in a good agreement with 
our analytical expectation equation 26. The small deviation from the analytical relation comes 
from our assumption that both the magnetic field and density will remain constant during the 
simulation. But the real situation is that both the density and magnetic field increase with 
time during the simulation. By taking the time evolution of both density and magnetic into 
account (inside the disk), the equation 25 becomes 

T e jection OC OC . (27) 

B z V A 

Both the density and B z component of the magnetic field increase steeply in case of strong 
initial magnetic field. On the other hand, the density almost remains constant in case of initial 
weak magnetic field while B z increases. Then as B z increases with time the ejection time 
becomes shorter than if we consider constant initial magnetic field. But the time evolution of 
the density will make an opposite effect in that case of strong initial magnetic field i. e. the 
effect of the increase of B z will be weakened by the effect of the increase in p. On the other 
hand, in the case of initial weak magnetic field, the effect of the increase of the density is too 
small to weaken the effect of the increase of magnetic field. Hence the evolution effect of B z is 
more prominent in case of weak initial magnetic field. So we think that the evolution of both 
density and magnetic field B z may explain the discrepancy between the analytical dependence 
of the ejection time and the initial magnetic field. Consequently, when we consider the average 
magnetic field instead of the initial magnetic field in Figure 11a, the agreement is more better 
as we can see from Figure lib. In that case the best fit is 

T r-i- Z? — 0-36 

average,eje ^mg 

5. Time averaged Velocities and Mass Outflow Rates of Jets As a Function of 
Magnetic Energy 

There is also an interesting question regarding magnetically driven jets: how the veloc- 
ities and mass outflow rates depend on the strength of the magnetic field. Figure 12a shows 
the time averaged jet velocities (Vj et ave ) as a function of the initial magnetic field E mg . Vj et av c 
is defined by 

Vjet ave = J q J q V z {r,z=z p )drdt j l$T orhi U ■ (28) 

where T e is the end of the simulation time, T^ut is the time of one disk rotation and 10T or bu is 
the sum of the time over which the speed of the jet is averaged. The figure shows that the jet 
velocity is 10 _1 smaller than that of order the Keplerian speed for a wide range of E mg in the 
disk and that its dependence on E mg as shown in Figure 12a is approximately. 

T/ rp 0.15 

''jet ave ^mg 
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Figure 12b shows the maximum jet velocities (Vj e t max) as a function of the initial mag- 
netic field E mg . The sold line shows 

V rsr F °- 17 

''jet max ^mg 

Figure 12c and 12d show the time averaged jet velocities (Vj e t avc ) and maximum jet 
velocities respectively as a function of E mg /M w . The sold line shows 

Vjet ave (^mg/ M w ) 

0.17 

Vjet max OC (E mg /M w ) ' . 

It is interesting that (in case of week magnetic field) this relation shows an approximate 
agreement with the scaling law of the Michel's solution (Michel 1969; Belcher & MacGregor 
1976); i.e., Voo = ($ 2 n F 2 /M w ) 1/3 , 

where is the poloidal velocity at infinity and $ = tq 2 Bq and VLp is proportional to 
the Keplerian angular velocity of the disk. 

Figure 13a shows the time averaged mass accretion rates of the disk as a function of the 
initial E mg . The figure shows that its dependence on E mg is about 



M a <xE mg 0S . 

Figure 13b shows the the time averaged mass ejection rates of our jets as a function of 
the initial E mg . Its dependence on E mg is approximately 

M w <xE mg 0A6 . 

Figure 13c shows the the time averaged toroidal magnetic energy as a function of the 
initial E mg . Its dependence on E mg is approximately 

E ♦ oc E a5 

Figure 13d shows the ratio of the time averaged mass outflow rate of the jet to the time 
averaged mass accretion rate of the disk as a function of the initial E mg . M w /M a is of order of 
0.1 or less 

M w /M a oc £ mg - a4 . 

5.1. Long time evolution of jet velocity 

Figure 14 shows the long time evolution of jet velocity which is calculated through the 
jet region r = — 1. Figure 14(a) shows the averaged jet velocity and Figure 14(b) shows the 
maximum jet velocity. Both the averaged jet velocity and maximum jet velocity never reach 
steady state and both of them have the same character of the variation. In case of strong initial 
magnetic field the jet velocity has the maximum value earlier than the weak initial magnetic 
filed case, but after the first 5 orbits the jet velocity becomes similar for all initial magnetic 
field. After the first 10 orbits the jet velocity for the initial weak magnetic field becomes higher 
than that in case of strong initial magnetic field case. 
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We notice that in all models the jet velocity decreases severely after the first three orbits. 
We think that the high velocity of the jet in the first stage of evolution is result from both the 
effect of initial and boundary condition. Also, at the first stage of the simulation the position 
of the ejection point of the jet is near the gravitational center. With the simulation going on 
the density and pressure of the central region of disk increases which leads to the decrease or 
stopping the accretion within that region. As a result for that the position of the ejection point 
of the jet becomes far from the gravitational center so that the jet velocity decreases. Also, The 
effect of softening gravitational potential near the gravitational center leads to the decrease of 
the jet velocity. While the simulation goings on, the magnetic field lines are trapped inside the 
softening gravitational potential. Consequently, they lose their angular momentum and their 
rotational velocities become very small which leads to the decrease of the jet velocity. 

5.2. Jet driving forces 

Fig. 15 shows the time evolution of powers of the jet i.e. Poynting flux F p j, kinetic flux 
Ffcj and enthalpy flux F en j which are described as: 



F j = - / 2irr—(E x B) z dr at z = 4, (29) 

JO 47T 

Fkj = / 2nrpv 2 v z dr at z = 4, (30) 

J 

F en ,j= [ 2iirpv z (—^—--)dr at z = 4, (31) 



,7-lp 

where z is z-coordinates for jets to pass through. 

The Poynting flux F p j plays the dominant effect in driving and accelerating the jet at 
the initial evolution of our simulation of a strong magnetic field case, (model 1) E mg = 2 x 10~ 4 
unit t — 50, after that it is suddenly becoming very weak. 

In case of initial weak magnetic field, model 5, E mg = 5 x 10 -6 , during the initial stage 
of the evolution, the enthalpy flux is the dominant one until t ~ 50, and both kinetic flux 
and Poynting flux have nearly the same value. After t ~ 50, the dominant energy flux is the 
Poynting flux and the other two fluxes have nearly the same value. At the last stage of the 
evolution the Poynting flux decreases again and the enthalpy flux becomes the dominant one. 
In case of very weak initial magnetic field , model 9, E mg = 2 x 10~ 7 , the enthalpy flux is the 
dominant one until the last stage of the simulation except after t ~ 225 it decreases sharply. 

Figure 16 shows the relation between the average enthalpy flux F en j ave , kinetic flux 
Fkj ave and Poynting flux F p j avc and the initial magnetic field strength. The average is defined 
as; 

F pJ ave = - J*" £ 27Tr ^( E x B) z drdt j 10T crit at z = 4, (32) 
Fkj avc = [ [ 2nrpv 2 v z drdt / 10T crit at z = 4, (33) 



o 
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F en ,j avc = J q J q 2nrpv z \^—^pj drdt J lOT crit at z = 4, (34) 

where T e is the end of the simulation time and 10T CT .j 4 is the sum of the time over which the 
flux of the jet is averaged. Figure 8 shows that the averaged enthalpy flux have the same 
value whatever is the initial magnetic field strength, whereas both the averaged kinetic flux 
and averaged Poynting flux increase with increasing the initial magnetic field. 

Kudoh & Shibata (1997a) showed that the dominant energy of a jet depends on the 
strength of the magnetic field. When the poloidal component of the magnetic field is B p oc r~ 2 , 
the fast magnetosonic point appears far from the Alfven point and the dominant energy of 
the jet is the Poynting flux. In our simulations, the initial magnetic field is uniform. In such 
models, the fast magnetosonic point locates far from the Alfven point (Kuwabara et al. 2000). 

6. Radial jet Structure 

Figures 17, 18 and 19 show the radial dependence of mass and energy flux of the jet. The 
mass flux definition is described by equation 18, the toroidal magnetic energy is described by 
equation 19 and the Poynting flux, kinetic flux and enthalpy flux are described by equations 32- 
34 respectively. Figure 17 show the radial profile of the density, poloidal velocity and toroidal 
field. The radial dependence of density shows that the peak density (which defines the jet) 
in models 1 and 4 (stronger initial magnetic field cases) is greater than that of model 6 (a 
weaker initial magnetic field case). The radial dependence of the poloidal velocity shows that 
the collimated flows show higher velocities closer to the disk axis. The radial dependence of 
the toroidal magnetic field shows that the maximum value is close to the disk axis as expected 
for collimated jet. Figure 18 shows the radial profile of the Poynting, enthalpy and kinetic flux. 
We plot the time averaged flux in annular rings around the symmetry axis (z-axis) to show 
the radial dependence of these energy flux. We show these radial dependence also by taking 
different initial magnetic field strength to show its effect. In case of model 4 the collimation of 
both mass flux and kinetic flux is clear. In figure 18 we calculated the flux at fixed height (z=4). 
Next we will calculate the flux at different heights. At each height we calculate the position of 
the maximum of each energy flux. Figure 19 shows the r-coordinate versus of z-coordinate of 
the maximum of mass flux, toroidal energy and kinetic flux for models 1, and 3. This figure 
shows how the maximum of each flux at different hight progresses in z-direction. From this 
figure we can notice that the collimation degree is different with the different kind of the flux 
and the initial magnetic field. In case of strong initial magnetic field the collimation is achieved 
earlier than weak case. Also the collimation is achieved after some height over the accretion 
disk. 
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7. Summary and discussion 



In this paper we have shown that long term 2. 5- dimensional MHD numerical simulation 
of magnetized accretion flows leads to an intermittent jet like outflow which never reach steady 
state. Our simulation is an extension of the works of Matusmoto et al. (1996) and Kudoh et al. 
(1998). Both of them performed time-dependent 2.5-dimensional MHD numerical simulations 
of jets from accretion disks including the dynamical processes within the disk. They showed 
that the ejection mechanism of the jets is the same as that in steady theory; i.e. the centrifugal 
force along the poloidal field line accelerates the jets within an Alfven radius and, above the 
Alfven radius, the jet is accelerated by magnetic pressure. The simulations for both of them 
last only for one inner disk rotation. What does happen for the characteristics of the jets if the 
simulations become long?. The answer of this question is given in this paper. 

7.1. The Physical Meaning of Time Evolution for M w , M a , and E mgt 

Cosmic jets are ubiquitous, being quite often associated with new-born stars, X-ray 
binaries, active galactic nuclei and gamma-ray bursts. In all such cases , jets and disks seem to 
be inter-related. Not only jets need disks in order to provide them with the ejected plasma and 
magnetic fields, but also disks need jets in order that the accreted plasma gets rid of its excess 
angular momentum to accrete. Observationally, there has already been accumulated enough 
evidence for such a correlation. For example, in star forming regions an apparent correlation is 
found between accretion diagnostics and outflow signatures (Hartigan et al. 1995). Hence, our 
current understanding is that jets are fed by the material of an accretion disk surrounding the 
central object. The mean controller between the mass accretion rate M a and the mass ejection 
rate M w is the magnetic field strength. Matusmoto et al. (1996) studied the dependence on 
the initial magnetic field strength. They showed that the ratio of the mass ejected as jet to the 
total mass of the flux tube was about 10% of the accreting mass. Figure 5 and figure 11 show 
that time evolution of the mass outflow and toroidal magnetic field outflow is quasi-periodic 
and the periodicity of the jet can be related to the time needed for the initial magnetic field 
to be twisted to generate toroidal field. Sato et al. (2003) found that these periodicities are 
around 2n. The toroidal magnetic field energy increases as magnetic field line is twisted and 
accumulated, because the magnetic field line is dragged with infalling and rotating gases. Then 
the piled up energy is released by magnetically driven outflow that is triggered by magnetic 
reconnection so that the mass outflow is intermittent and has some periodicities. 

Figure 13d show that the average M w and M a (averaged over 10 rotations) is closely 
related. The magnetic field strength is the controller of the relation between M w and M a . 
When the initial magnetic energy is weak the ratio between M w and M a is high. The ratio 
decreasing with increasing the initial magnetic field. Pelletier & Pudritz (1992) stated that if 
the disk magnetic field is reduced, there is a very large increase in the mass-loss rate in the wind. 
The point is that since a slower wind is driven in weaker disk fields, one must provide much 
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more mass in order to carry off the same amount of disk angular momentum. At strong initial 
magnetic energy field case the ratio reach a constant value. These dependences are consistent 
with the results of the steady solution (Kudoh & Shibata 1997a). 

7.2. The Dependences on the Initial Magnetic Field Strength in the Weak field Case 

Kudoh & Shibata (1997a) derived the dependence of the mass outflow rate (M w ) and 
mass accretion rate (M a ) on the initial magnetic field strength using a semi- analytical method, 

M w oc E m V 2 , (35) 

M a oc£ mg . (36) 

Furthermore, Michel's scaling law (Michel 1969) is written as 



On the other hand, we obtain 

M w ,a,ve rp -0.4 /oq\ 
-'"o,ave 

V jet> a Ve oc E^ 5 . (39) 

Our results shown in Figure 12 and Figure 13 are consistent with not only the steady 
solutions but also the results of both Kudoh et al. (1998), Kato et al. (2002) and Kigure & 
Shibata (2005) However there is a remarkable difference from the results of Kudoh et al. (1998), 
Kato et al. (2002) and Kigure & Shibata (2005). They show only maximum values for M w , 
M a and V z , while we examined time-averaged values. We find that V z ave and M w are 0.1 times 
smaller than V z max and M w max in the case of Kudoh et al. (1998). 

8. Conclusions 

We have performed time- dependent, 2.5-dimensional, axially symmetric, MHD numerical 
simulations of jets for many orbital periods. Our simulations solve the dynamical process of the 
accretion disk for longer time case than that of the previous simulations. Because we wanted 
to know whether the jet continued to be ejected and become steady (or quasi-steady) state, we 
investigated time variations of M w , M a , and E mgt . We also wanted to know the dependence of 
M w ave , M a ave , and V z ave on the initial magnetic field strength and compared our results with 
the steady theory and cases of Kudoh et al. (1998) and Kato et al. (2002). We summarize our 
results as follows: 

1. In all models, the ejection of jets is intermittent. 

2. The ejection of jets has a period which related to the toroidal field formation. 
There are also relationships between M w and M a , and between M w and E mgt ; specifically, 

their time variations are similar. 
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3. The dependence of M w , M a , and E mgt on the strength of the initial magnetic field 
is consistent with those in the steady theory and cases of Kudoh et al. (1998) and Kato et al. 
(2002). In all cases, however, V z ave and M w are 0.1 times smaller than V z max and M w max in 
the case of Kudoh et al. (1998) and Kato et al. (2002) and Kigure & Shibata (2005) 
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TABLE 1 



UNITS FOR NORMALIZATION 





Normalization 


Physical Quantities 


Unit 


t(time) 




r,z (length) 


ro 


^(density) 


Po 


p(pressurc) 


PoV£ 


v (velocity) 


Vko 


B (magnetic field) 




e (specific internal energy) 


Vio 



Table 1. The unit length ro = (L^/GM )'l/(l — 2a) is the radius of the density maximum in the initial disk. The 
unit velocity Vfc = (GM/ro)(l/2) is the Keplerian velocity at (r, z) = (ro,0). The unit density po is the ini- 
tial density at (r, z) = (ro,0). It is assumed that a = 0.45, in our study (see eq.12 for the definition of a) 

TABLE 2 
MODEL PARAMETERS 



model E mg E t h time 

1 2 x 1CT 4 0.018 80.0tt 

2 8xl0" 5 0.018 22.1tt 

3 5xl0" 5 0.018 31.7tt 

4 2xl0" 5 0.018 31.2tt 

5 5 x 10~ 6 0.018 80.0tt 

6 2 x 10~ 6 0.018 37.7tt 

7 8xl0" 7 0.018 41.3tt 

8 5 x 10~ 7 0.018 34.8tt 

9 2 x 10~ 7 0.018 80.0tt 
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Fig. 1. Temporal evolution of an ideal accretion disk threaded by a poloidal magnetic field. Color levels 
represent density level (upper panel) and temperature level (bottom panels) while solid lines stand for 
poloidal magnetic filed lines. The time unit labeling each snapshot is for t ~ 2n ~ 6.28 corresponds to 
one Kcplerian orbit at (r, z) — (1,0). After a few rotations, outflows are escaping from the disk and the 
outflow and accretion remain until the last moment of our simulation time. The square refers to the region 
analyzed by Kudoh et al 1998 
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Fig. 2. Temporal evolution of an ideal accretion disk threaded by a poloidal magnetic field. Color levels 
represent density level (upper panel) and temperature level (lower panel) while solid lines stand for poloidal 
magnetic filed lines. The time unit labeling each snapshot is for t ~ 2n ~ 6.28 corresponds to one Keplcrian 
orbit at (r,z) = (1,0). After a few rotations, outflows are escaping from the disk and the outflow and 
accretion remain until the last moment of our simulation time 
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Fig. 3. Temporal evolution of an ideal accretion disk threaded by a poloidal magnetic field. Color levels 
represent density level (upper panel) and temperature level (bottom panels) while solid lines stand for 
poloidal magnetic filed lines. The time unit labeling each snapshot is for t ~ 2ir ~ 6.28 corresponds to one 
Kcplerian orbit at (r,z) = (1,0). After a few rotations, outflows are escaping from the disk and the outflow 
and accretion remain until the last moment of our simulation time 



24 



f ig4 . eps 



Fig. 4. Ratio of B v to B p along r = 0.225. at different three time shot, for different magnetic field energy. 
The twisting magnetic field become more prominent in case of weak magnetic field. 
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Fig. 5. Time variation of the mass outflow rate for three models, the out flow is still intermittent after 
long time simulation for different initial magnetic field energy 
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Fig. 6. Time variation of the mass accretion rate for three models, the accretion is still intermit- 
tent after long time simulation for different initial magnetic field energy. One Keplerian orbit at 
(r,z) = (1,0) ~2tt~ 6.28 time unit. 
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Fig. 7. Accretion rate, averaged between r = 0.52 and r = 10, against time for different initial magnetic 
field energy. One Keplerian orbit at (r,z) = (1,0) ~ 2ir ~ 6.28 time unit 
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Fig. 8. Accretion rate as a function of radius, for different initial magnetic field, averaged over 10 orbites 
for models 1, 3, 4 and 6 
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Fig. 9. Time variation of the toroidal energy for three models, the outflow is still intermittent, similar to 
mass out flow, after long time simulation for different initial magnetic field energy. One Keplerian orbit 
at (r,z) = (1,0) ~ 2tt~ 6.28 time unit 
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Fig. 10. The mass flux for four models as a function of time. The numbers correspond to the maximum 
flux at which we measure the corresponding time (peak time ) and plot this time with the initial magnetic 
held in fig. 11 



31 



f igllpp . eps 



Fig. 11. In this figure we show the relation between the average time interval between the mass ejection 
peaks and the initial magnetic field strength (a) and the average magnetic field (b). Solid line corresponds 
to the best fit to the numerical result and the dotted line corresponds to the analytical relation ( equation 
26). 
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Fig. 12. (a) Time averaged velocities (Vj et a ve) as a function of E mg .(b) V z max as a function of E mg , (c) 
Time averaged velocities(Vj- e t avc) as a function of E mg /M w . (d) Vj e t max as a function of E mg /M w . 
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Fig. 13. (a) Time averaged mass accretion (M a ) rate of the disk as a function of E mg . (b) Time averaged 
mass ejection (M w ) rate of the jet as a function of E mg . (c) Average toroidal ejection (E mgt ) as a function 
of £mg- (d) The ratio of the time averaged mass outflow rate of the jet to the time averaged mass accretion 
rate of the disk (M w /M a ) as a function of E mg . 
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Fig. 14. The jet velocity for three models as a function in time, (a) The averaged jet velocity (averaged be- 
tween r = and r = 1 ) at height z = 4 as a function in time, (b) The maximum jet velocity at z = 4 as a func- 
tion in time (averaged between r — and r = 1. One Kcplcrian orbit at (r,z) = (1,0) ~ 27r~ 6.28 time unit 
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Fig. 15. Time evolution of powers of the jet; Poynting flux, F p j (solid line), thermal flux F en j (dashed 
line), and kinetic flux Fkj (dotted line) for different initial magnetic field 
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Fig. 16. The time average value of enthalpy, kinetic and Poynting flux as a relation with the initial 
magnetic field strength. 
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Fig. 17. Radial profiles of average physical quantities. The cuts were taken at z = 4. The radial de- 
pendence of density for each initial magnetic field shows that the peak density (which defines the jet) in 
the first two initial magnetic field is greater than the last weak initial magnetic field case. The radial 
dependence of the poloidal velocity shows that the collimated flows show higher velocities closer to the 
disk axis. The radial dependence of the toroidal magnetic field shows that the maximum value is close to 
the disk axis as expected for collimated jet. 
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Fig. 18. This figure shows the flux outflow of (a) mass flux ,(b) kinetic flux,(c) toroidal magnetic energy, 
(d) enthalpy flux, and (e) Poynting flux at z=4 as a function of r-coordinate. The collimation becomes 
more clear in case of model 4 in both mass flux and kinetic flux 
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Fig. 19. The lines represent the progress of the maximum of different fluxes in z-direction. The collimation 
begins at low-z in case of strong initial magnetic field. 
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